clear; close all; 
addpath('../utilities');

rng('default')

% characteristics of 2008 disaster
size_2008   = -0.0214; % mean 2008 size (simple net growth rate)
margin_2008 = 0.001;


% load consumption data
dc          = xlsread('../data/consumption.xlsx','H15:H104');
start       = 19; % 1947-48 growth rate
data_mom    = [mean(dc(start:end)),std(dc(start:end)),skewness(dc(start:end)),kurtosis(dc(start:end)),corr(dc(start:end-1),dc(start+1:end))]';


p       = 0.878776;

mu1     = 0.003060;
mu2     = 0.000453;

sig1    = 0.001500;
sig2    = 0.002768;
sig3    = 0.007034;

DD           = 44; % mean disaster duration in month
i            = 1:100000;
duration_fun = @(p)sum((1-p)*p.^(i-1).*i);
p3          = fzero(@(p)duration_fun(p)-DD,0.99);

ED          = 0.834185141; % mean drop in consumption (simple growth rate)
i           = 1:100000;
mean_fun    = @(m3)sum((1-p3)*p3.^(i-1).*exp(i*m3+i*sig3^2/2));
mu3         = fzero(@(m)mean_fun(m)-ED,-0.0038);

pd          = fzero(@(x)solve_P(x,x,p,p,p3)-(44/12)/(2019-1929),1e-3); % This is pi in the paper.

mu_C        = [mu1;mu2;mu3];
sig_C       = [sig1;sig2;sig3];

P = [ p*(1-pd)      (1-p)*(1-pd)  pd;
      (1-p)*(1-pd)  p*(1-pd)      pd;
      (1-p3)/2      (1-p3)/2      p3];
  
  

% generate '2008 recession' paths
Dsize        = size_2008; % mean drop in consumption (simple growth rate)
Dmin         = Dsize-margin_2008; % range of cummulative drop
Dmax         = Dsize+margin_2008;

Tsim         = (2019-1995)*12 - 6; % Jan/1996 -- Jun/2019
Psim         = 7000;
numSamples   = 5000;

states       = NaN(Tsim,Psim);
epsC         = NaN(Tsim,Psim);
delC         = NaN(Tsim,Psim);
max_drop     = NaN(Psim,1);

for i=1:Psim
    
    found = false;
    for pp=1:1000
        
        if i==1
            s0 = 1;
        else
            s0 = states(end,i-1);
        end
        
        [~,states_i] = sim_markov_chain_fast(P,Tsim,1,s0,1:3);
        epsC_i       = randn(Tsim,1);
        delC_i       = mu_C(states_i) + sig_C(states_i).*epsC_i;
        
        max_loss = 0;
        cum_loss = 0;
        for t=1:Tsim
            if delC_i(t)<0
                cum_loss = cum_loss + delC_i(t);
            else
                if cum_loss<max_loss
                    max_loss = cum_loss;
                end
                cum_loss = 0;
            end
        end
        
        if (exp(max_loss)-1)>=Dmin && (exp(max_loss)-1)<=Dmax
            states(:,i) = states_i;
            epsC(:,i) = epsC_i;
            delC(:,i) = delC_i;
            max_drop(i) = max_loss;
            found = true;
            break
        end
        
    end
    if ~found
       break
    end  
end

% do not start with a disaster, so that we can initialize the firm distribution from a post-war (no disaster) simulation
ix = states(1,:)~=3;
states = states(:,ix);
epsC = epsC(:,ix);
delC = delC(:,ix);
states = states(:,1:numSamples);
epsC = epsC(:,1:numSamples);
delC = delC(:,1:numSamples);

save aggregate_path_samples_long states epsC delC







